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THE ONE-DIMENSIONAL TIME-DEPENDENT INTERACTION OF A SATELLITE 
WITH THE IONOSPHERIC PLASMA 


SUMMARY 


The time -depen dent equations describing the 
one-dimensional interaction of a streaming plasma 
with an infinite flat plate were derived and solved 
numerically by a modified Euler method. Solutions 
were obtained for: (i) the case of a plate with no 
initial surface charge at an altitude of 200 km and 
(2) the transition from equilibrium conditions at 
200 km to an equilibrium state at 3000 km. 

In general the surface charge density was found 
to overshoot to a peak value early in the interaction, 
as a result of the initial rapid collection of electrons, 
and slowly decay toward an equilibrium value as the 
effect of the ions became more pronounced. At 
200 km the sheath thickness was found to be 0. 33 cm 
and the relaxation time 1.42 ^ts, while at 3000 km 
the sheath thickness and relaxation time were 7. 3 
cm and 30. 5 jus respectively. 


INTRODUCTION 

Purpose 

The net result of the interaction of a body, 
whether moving or stationary, and a plasma, in 
which the body is immersed, is a charge accumu- 
lation on the surface of the body. The electric 
potential associated with this charge in turn affects 
the surrounding plasma, creating a space charge in 
the disturbed or sheath region near the body. 

In this report, a theoretical treatment of the 
mechanisms that produce the electric potential 
associated with a body moving through a collision- 
less ionospheric plasma is undertaken. Special 
attention is given to the time dependence of the 
electric field with the intention of investigating the 
behavior of the various sheath parameters as the 
body moves rapidly through non uniform regions of 
the plasma. A system of equations that describes 
a one-dimensional model of the problem is established, 
and numerical methods are invoked to find a solution. 


Background 

The realm of plasma physics dealing with 
satellite-ionospheric interactions has grown in 
interest along with the ever increasing application 
for and complexity of earth satellites. The sheath 
effects mentioned above can significantly alter the 
performance of onboard instrumentation and, 
therefore, must be taken into account in the planning 
and interpretation of sate llite -borne experiments. 
Correspondingly, a great volume of work has been 
done in this area. 

To briefly review development in the field, we 
should first consider an early paper by Jastrow 
and Pearse in which the effects caused by the 
motion of a body in a plasma are studied, assuming 
a spherically symmetrical potential distribution 
and a uniform, undisturbed ion density LlJ. These 
assumptions have since been found to be incorrect 
so that the above approach is good only as a first 
approximation. In a later paper, Kraus and Watson 
used kinetic equations to calculate the potential 
and ion number density around a point-like charge 
moving through the ionosphere [2J. This would 
imply that the body radius is much smaller than 
the Debye radius. However, it has been found that 
the Debye radius in the ionosphere is considerably 
smaller than the dimensions of a typical satellite. 

A third paper, published by Davis and Harris and 
dealing with the same problem, is of interest, 
because numerical methods were used to solve 
simultaneously the equations of potential and the 
motion of ions near a body moving in a plasma [3J . 

Although these papers provide a good foundation 
in the area of satellite-ionospheric interactions, 
none of them attack the time -dependent problem. 

In fact, in searching the literature for information 
on the problem treated herein, it was found that 
relatively little had been done toward the develop- 
ment of a suitable time -depen dent theory of the 
satellite-ionospheric interactions, most of the 
work having dealt with all aspects of the steady- 
state interaction for various vehicle geometries. 

Some work of a time -dependent nature has been 
done by G. S. Kino and associates at Stanford 



University [4, 5, 6, 7] . However, this work deals 
with the problem of transmitting rf waves through 
a discharge plasma and does not apply directly to 
the problem treated herein. Getmantsev and Denisov 
have also done work of a time- dependent nature in 
which they found a plasma disturbance caused by 
high-frequency electromagnetic field fluctuations 
near cylindrical bodies slowly moving through a 
plasma f8]. This work applies only to slowly 
moving bodies. This is not a good assumption for 
satellites and, further, would change the nature of 
the ion interaction since their thermal energy 
would become important. 

Two particular cases of the time -depen dent 
satellite-ionospheric interaction will be treated in 
detail in this report. First, there is the case in 
which the satellite has no initial surface charge and 
is assumed to be at zero potential. Assuming an 
equal density of ions and electrons in the undisturbed 
plasma, the more mobile electrons will initially 
collect on the surface of the satellite at a greater 
rate than the ions. This will create a negative 
potential on the satellite that will, in turn, slow the 
incoming electrons and accelerate the ions, thereby 
increasing their relative rate of incidence. Hence, 
as the ion flux at the surface builds, the potential 
will become more negative at a decreasing rate until 
an equilibrium condition is reached in which the ion 
and electron fluxes at the surface of the satellite 
are equal. 

In the second case treated herein the interaction 
is initially taken to be at some preestablished 
equilibrium state. A step change is then made in 
the plasma parameters and the interaction moves 
to a new equilibrium state. 

The general behavior of the interaction in this 
case is similar to that described previously with 
the exception that there is an initial surface charge 
on the satellite and, hence, an electric field that 
acts on the plasma. The effect of the electric 
field on the plasma produces the equilibrium state 
of the previous case. From this point, the electric 
field increases or decreases to become compatible 
with the new parameters. 

Approach 

In order to attack the problem, it was necessary 
to make several assumptions. First, by talcing the 
body or satellite under consideration to be an infinite 
flat plate, oriented perpendicularly to the plasma 


flow, the problem is effectively reduced to one 
dimension with a corresponding reduction in com- 
plexity of the equations. This can be interpreted 
physically to mean that only a very small region at 
the stagnation point on the satellite will be con- 
sidered. Consequently, the electric field, and 
other sheath parameters, can be observed only 
in the forward direction . However, the problem 
is still physically meaningful, since it is conceivable 
that an infinite flat plate would be a reasonable 
approximation of a small probe located at the 
stagnation point of a satellite. 

To facilitate the formulation of the problem, we 
take the inverse of the actual physical situation and 
let the plasma stream into a stationary plate with 
the satellite velocity. This will not affect the solu- 
tion since it is the relative velocity of the satellite, 
or plate, and the plasma that are important. It 
should also be noted that any experimental verifica- 
tion in earthbound laboratories will require this 
arrangement. Since the ion mean thermal velocity 
in the ionosphere is an order of magnitude less than 
typical satellite velocities, the random thermal 
motion of the ions will be neglected. The ions, 
then, will all have the same velocity at any time and 
position. The mean thermal velocity of the 
ionospheric electrons is much greater than typical 
satellite velocities, and we therefore assume them 
to have a Maxwellian velocity distribution. Con- 
sidering the small ratio of electron to ion masses, 
it is reasonable to assume an instantaneous relaxa- 
tion of the electrons. 

In this treatment we omit the effects of the 
earth* s magnetic field, neutral particle collisions, 
and corpuscular radiation (photoelectric effect) . 

Since we are dealing with a collisionless plasma 
in a very small region near the frontal surface of 
a satellite, these components of the interaction will 
produce relatively small effects compared with the 
electric potential resulting from differential charge 
flux to the surface. 

Two further facts about the ionospheric plasma 
should be brought out before going further into the 
treatment of the problem. First, as a result of 
neglecting the ion thermal velocity, we effectively 
assume zero ion temperature, and therefore only 
the electrons will be affected by temperature. Any 
further references to temperature will, therefore, 
be made to the thermal energy of the electrons. 
Secondly, the undisturbed, ionospheric plasma is 
maeroscopically neutral so that at any point far 
from the plate the electron and ion number densities 
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will be equal. Some feeling for the validity of these In terms of electrical charge, the ionosphere 

and the other assumptions made in this chapter can can be considered to be a two-component gas con- 
be obtained by observing the ionospheric parameters. sisting of singly charged ions and electrons. A 
Some typical values are given in Table i. separate Boltzmann equation must be applied to each 

TABLE 1. IONOSPHERIC PARAMETERS 




Altitude (km) 


Parameters 

200 

300 

3000 

Satellite Velocity (km/ sec) 

7.78 

7.73 

6 . 52 

Ion-Electron Density (m“ 3 ) [9] 

(3-50) x IO 10 

( 10 - 20 ) x 10 11 

7 x 10 9 

Temperature (K° ) [9] 

450-800 

1000 

4000 

Ion-Electron Mean Free Path (m) [ 9 J 

90 

70 

3 x 10 4 

Neutral Mean Free Path (m) [9] 

80 

1000 

2 x 10 12 

Debye Length (cm) [9] 

0 . 2-1 

0. 14-0.7 

4 

Ion Composition [ 10 J 

NO + , 0 2 \ 0 + 

O 4 (98^), N 4 

+ 

H 

e 

Average Ion Mass 

24 

16 

4 

V /V. 
s 1 

13 

11 

4 


DERIVATION OF EQUATIONS 
Time- Dependent Equations 

The ionospheric plasma must be treated as a 
collection of discrete particles rather than as a 
continuous medium. It follows, then, that fluid 
flow equations used in hydrodynamics and aero- 
dynamics cannot be applied to this type of problem. 
Nor is it feasible to follow dynamic trajectories of 
individual particles, except in very limited cases. 
These difficulties can be circumvented, however, by 
use of the Boltzmann equation, which provides 
accurate statistical information about the distribution 
function and average expected values of quantities 
describing particle behavior [11]. In its collision- 
less form, the Boltzmann equation is 


species of a multicomponent gas [12]. We will 
therefore have two equations, one that describes 
ion velocity distribution and another for electron 
velocity distribution. 

As pointed out previously, the electron thermal 
velocity is much greater than the satellite velocity, 
and it is therefore reasonable to consider this 
component of the ionospheric plasma to be at 
equilibrium. In this case, the electron behavior 
can be described by the Maxwell-Boltzmann distri- 
bution function [ 13 } 14 J . i n one dimension, this 
function is of the form, 


/ m 1/2 

me 2 

rf e \ 

e $(x,t) 

\ 27rkT ) 6Xp 

” 2kT “ kT 

N e' 

e e 


( 2 ) 


af af a 

at ax ac 


(fF) = 0 


( 1 ) 


where f is the distribution function, c is velocity, 
and F is acceleration resulting from external 
forces. 


where N is the total number of particles in the 
system, k is the Boltzmann constant, and <f> is 
the potential energy of the system. This is an exact 
solution to the steady-state Boltzmann equation. Its 
use as a time-dependent solution assumes instan- 
taneous relaxation of the electrons. 
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We now consider the ions. As stated previously, 
the ion thermal velocity is much smaller than the 
satellite velocity and can be neglected. The ions 
will then have a uniform velocity distribution in 
space and time, which can be written in the form, 

f.(x.t,c) = n.(x,t) 6 [ c + u.(x,t) ] (3) 


where u. is the average ion velocity and depends 
on x and t, and n. is the ion number density 
and also depends on x and t. 

Multiplying the Boltzmann equation (1) by a 
constant or any function of ion velocity and inte- 
grating over velocity space, we can obtain the 
equations of transfer, which are more easily 
solved than the Boltzmann equation [15]. Using the 
ion mass, im, as the multiplier and performing 

the appropriate integration, we obtain 


0 


f 


m. 

l 


df. df. 

-L + c — 
dt dx 


d 

dc 


(Ff.) 


dc 


0 


Notice that the third term reduces to 
Fn.6(c-uJ (° , which is zero at both limits. 

l i 

The remaining terms have values only when 
c = -u.. The foregoing expression therefore 

reduces to 


a 

at 


(n i ) + ir 



= 0 


or 



d 

dx 




n.(x,t) u. (x,ti 


= 0 . 


(4) 


This is the conservation of mass, or continuity 
equation for ions. 


By multiplying the Boltzmann equation (i) 
by 1/ 2 m.c 2 and integrating over velocity space 

in a similar manner, we can obtain the conserva- 
tion of energy equation for ions. The integration 
is carried out as follows: 


or 


0 af. 0 af. 0 

f dc + f c dc + f — - (Ff.) dc = 0 . 
J 9t J dx j ac i 


Since x and t are independent variables, the 
partial derivatives of these variables can be 
brought outside the integral. Performing thi^ 
operation and replacing f by expression (3). we 
obtain 


0 


/ 


m.c 2 

i 


af. af. 

i i 

" + c — — 

at ax 



(Ff J 


dc = 0 



0 

f c 2 6(c + u.) dc 

J l 

— oo 



0 

n. J c 3 6(c - u.) d< 

— C 50 


d_ 

at 


o 

n. J 6(c + u.) dc 

— OO 



_a_ 

ac 


6(c + u.jj 


dc 


0 


+ 


_a_ 

ax 


o 

n. J c6(c + u.) dc 

— ' oo 


+ Fn. 


..j —[fitc + v] 


dc = 0 


The acceleration, F, results from the electric 
field associated with the plate and is therefore 
qE/nr , where E is a function of x and t. The 

acceleration, F, can therefore be taken out of 
the integral since in this case it is not a function of 
velocity. This yields 


4 



a 

at 


[n.(-u.) 2 +-~- n. (-u.) 3 

1' 1 ax 1 1 


+ Fm | c 2 6(c + 


- J 2c 6 (c + m) dc > = 0 


Note that the term c 2 6(c - u.) |° is zero at both 

1 -CO 

limits, and the integral is 2(-u,.). The preceding 
expression now simplifies to the form, 

Tr[ n i (x * t) “h*- 4 *] -~k [ n i (x>t) u l (x>t) ] 

+ — E (x, t) n.(x,t) u. (x,t) = 0 
m. i ’ i 

i 


where we have replaced F by qE/nm . This 

equation can be simplified if we first expand the 
terms to obtain 


2n. u. (u.) + u 2 ~~~ (n.) - n.u. (u 2 ) 

i i at i i at x i i i ax ' i 


- u 2 •— — (n.u.) + — — En. u. = 0 
l ax i i m. ii 
i 


Rearranging terms we get 


2n.u. ~~ (u.) + \i 2 . (n.) - — — (n. u.)l 
i i at ' t i [ at v i ax i i 

- 2n. u? — — (u.) + En. u. - 0 

i i ax l m. i i 
i 


Now the term in brackets is exactly equation (4), 
which is equal to zero. We can also divide the 
remaining terms by 2m u to obtain 


d 

at 

U. (X, t) 

- u i (x ' t,_ ^r 

u. (x, t) 




_ 


+ — E(x,t) = 0 
m. 

l 


(5) 


which is the conservation of energy equation 
for ions. 

In the above derivations, the limits of integra- 
tion are taken from -°° to 0, because the ions 
move in the negative x-direction when traveling 
toward the plate. The 0 to +«> range is not included 
because all ions that strike the conducting plate are 
absorbed or neutralized. Since all ions initially 
stream toward the plate and none are reflected, 
there can be no positive component of the ion 
velocity distribution (Fig. 1). 

The two foregoing conservation equations 
contain the three dependent variables, m, u., and 

E. A third equation is, therefore, required to 
obtain a solution. This is provided by Gauss 1 
equation, which in its general form is 


V • E(x,t) - p(x, t)/c 0 . (6) 


In the above equation p is the net charge density 
per unit volume. Note that we have not used the 
wave equation because the potential retardation 
term, i/c 2 (9 2 <£/9t 2 ) , is assumed to be negligible 
for this case. 

The net charge density is obtained by adding 
the total charge density of ions and electrons as 
follows: 


P(x, t) 


q £ Zn. (x,t) - n^ (x, t) J 


(7) 


Here, Z represents the number of charges per 

ion. However, since we assume only singly 

ionized particles to exist, Z will have a value of 

unity. The number densities, n. and n , are 

l e’ 

average values. The ion average number density, 
n., must be obtained from the solution of the 

l 

equations; however, the electron average number 
density can be found from expression (2) as follows 
[16]: 
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CONDUCTING FLAT PLATE 


Figure 1. Plasma stream and plate configuration. 


OO 



The potential energy for electrons, <£, is equal 
to -q <p where cp is the electric potential. This 
can, in turn, be written in terms of the electric 
field as follows: 

E = - V(p = - dcp/dxi 


OO 0 

J E dx = - J d <p' = <fi 

X ( j ) 


N 


l 2?rkT e 


1/2 


-VkT e 

e 


2?rkT 

e 


m 

e 



-^ /kT e 

n = N e 
e 


Then the electron number density can be written 
in the form 


OO 

i^r f E(x>t) 

n e (x,t) = Ne e X ' 


(8) 


6 


Substituting equation (8) and equation (7) into 
equation (6), we obtain Poisson r s equation in the 
form 


t <E (U. i) 
Ot 


ii 

€ 


n. (0, t) ul (0, t) - N u e 


kT 

e 


CO 

/ 


E (x, t) cL\ 


9E(x,t) = JL 
9x e 0 


q 

kT 


n.(x, t) - N e 


OO 


J E (x, t) dx 

x 1 


(9) 


( 10 ) 

which is the desired boundary condition. 

The set of initial conditions will depend on the 
particular case being solved and will therefore be 
given later in the report rather than here. 


The boundary conditions for the system of 
equations (4), (5), and (9) can be described in 
terms of the electric field, which must go to zero 
for any set of initial conditions and time as the 
distance from the plate approaches infinity. At 
the plate, the value of the electric field will depend 
upon the initial conditions. However, its change 
with time can always be represented by the differ- 
ence of the ion and electron fluxes to the plate. 
Therefore, 


da 

dt 


q(n. u. - n u ) 
i l e e 


x 


0 


Nondimensional Time-Dependent Equations 

The set of equations derived in the previous 
section will be more manageable in the numerical 
analysis if they are nondimensionalized. We 
define nondimensional density, D = n./N; velocity, 

V = u./V ; distance from the plate, y = x/L; 

1 S A 

time, T = t / t, and electric field E = E/£, where 
N is the undisturbed ion number density, V g is 

the satellite velocity, and L, r, and £ must be 
determined. In terms of these variables, equation 
(4) takes the form 


The average values, n^ and u^, can be 
represented as follows: 


n 

e 



I Edx 

0 


This was shown previously in equation (8). 
Similarly, 


u 


e 



f dc 
e 


u e 


q 

kT 

e 


/ Ed x 
0 


where 

u = (kT /2irm j 1 / 2 
e e 


Now, from Gauss' Law it can be shown that 
E = o/€ 0 [17]. Applying this and the two 
expressions for n and u found previously, 
we obtain 


9D 

9T 


r V 

s 


L 



0 


Now if we let r = L/V , then the previous equa- 
tion becomes 


d_D(y,T ) T) 9V(y,T) 

9t ' 9y 


- V(y,T) 


9D(y,T) 

ay 


0 


( 11 ) 


which is the nondimensional continuity equation. 

Now consider equation (5). Upon substitution 
of the nondimensional variables, it takes the form 

V V 2 

_S_ _9V __ s. 9V +JL A 
r 9 T L 9y itl 

Now recall that t has been defined as L/V . 

s 

Making this substitution leads to the form, 


_ V ^V + 

9T 9y 



0 
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By defining a = (q£L)/ (mW^ 2 ) where a deter- 

mines the amount of effect the electric field has on 
ion velocity, we obtain the following nondimensional 
form of the energy conservation equation: 


my^L) M +a i (y>x) = o 

dT ay 


(12) 


Similarly, equation (9) takes the following form: 


9E 

ay 


q NL 
£ € o 


D - e 


kT f 
e y ! 


r a 

f E d y 


A 2/E(y,T)dy' 

8E( 8 y } = D(y, T) V(y,T) - V e ° 

(14) 

which is the nondimensional form of the boundary 
condition at the plate. 

The factors L, r, and 4 have all been 
defined. Now each can be expressed in terms of 
the plasma parameters. We already have 

L = n/ e 0 kT e /q 2 N. Using this expression for L 

in the expression for r, we obtain 

t = \/ e 0 kT e /q 2 NV g 2 . Similarly, 4 = ^NkT^/^o 

and a = kT /m.V 2 . 

e i s 


The exponent must be dimensionless, so we define 
4 - kT^/qL. Substituting this definition of 4 into 

the coefficient of the term on the right side of the 
equation yields q 2 NL 2 /(q)kT e ) , This can be made 

dimensionless by defining L = 4 i^kT^7 q*n. The 

above equation now becomes 


Steady-State Simplification 

The equations describing the steady-state 
interaction can be derived directly from the time- 
dependent equations (11), (12), (13), and (14) 
be setting all derivatives with respect to time 
equal to zero. Thus equation (11) becomes 


9E (y, T) 

ay 


/ E (y, T) dy 

D(y,T) - e y , (13) 


= 0 . (15) 

By dropping the time derivative in equation (12) 
we get 


d_ 

dy 


D(y,T) V(y,T) 


which is the nondimensional Poisson f s equation. 

a dV 

a E = V “ — 

Finally, equation (10) must be nondimen- dy 

sionalized. Making the variable substitutions, it 

takes the following form: Integrating both sides of this equation with 

respect to y we obtain 


A 

8E 

dT 


qr N V 

s 

€ 0 € 


DV - u/V g 





or 


Recall now that t - L./V so that the coefficient 

s 

of the term on the right side of the equation becomes 
(qNL/e 0 )/4, which is unity since 4 = kT^/qL and 

L = \l €„ kT e /q z N. Now if we let V = u/V g) the 

above equation becomes 


P(y) 


~ [ V 2 (y) - 1 
Za 


(16) 


where P denotes the nondimensional electric 
potential. 
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Equation (13) takes the form 



or in terms of the nondimensional, electric 
potential, P, it is 


d 2 P(y) 

dy 2 


e p (y) 


- D(y) 


(17) 


Equations (15), (16), and (17) can now be 
combined into a single equation for the electric 
potential. Notice that equation (15j implies that 
DV = constant. This constant is required to be 
unity by the boundary conditions on ion density 
and velocity at infinity. Then we obtain the solution 


D(y) - 1 /V (y) . (13) 

Equation (16) can be solved for V(y) to obtain 

\j 1-3 a P. Replacing V in equation (18) with 
this expression we obtain 


D - 1 M 1-2 cv P 


This expression can, in turn, be used to replace 
D in equation (17) to obtain 


The integral of the electric field, in this equation, 
is just the electric potential at the origin, or plate, 
which will be designated as P 0 . Making this 
substitution and solving for P 0 we find that 

P 0 = In (DV/V) , (20) 

which is the potential on the plate for a given set 
of plasma parameters. 

The solution of equations (16), (17), (18), 
(19), and (20) by numerical methods is simpler 
and less subject to error than a numerical solution 
of the time- dependent equations of the previous 
section. Their solution will therefore provide a 
check on the degree of approach to the steady-state 
solution of the time -dependent equations. 


APPLICATION OF NUMER ICAL METHODS 
Modified Euler Method 

Given a differential equation and its solution 
at some initial or boundary point, the solution can 
be extrapolated to a neighboring point by use of the 
Taylor series. This is an infinite series of the 
form 


Y , - Y + M h + - 7 - M' h 2 + -r M” h 3 + . . . 
n+l n n 2 n 6 n 

( 21 ) 


<i 2 p (y) 

dy 2 


e P(y) - - 

\j 1-2 « P (y) 


(19) 


which provides an equation of P as a function of 
y and the plasma parameters. 

The boundary condition at the plate given by 
equation (14) can also be used in the steady-state 

A 

solution. Equating 8 E/at to zero, we obtain 


OO 

_ */ 
DV -V e ° 


A 

E dy 


0 


where M - (dY/dx) , JM' - (d 2 Y/dx 2 ) , . . . , all 
n n n n 

of which are evaluated at the point x - x 0 + nh, 

where h - x , - x . Although this method is 
n -1 n 

accurate if h is sufficiently small, it is not 
frequently used in this form since it usually requires 
more work than other methods. It is useful, 
however, in finding the first few points of a solution, 
which, as will be seen, cannot always be obtained 
from other methods. 

The method used herein is a modification of the 
Euler method. To describe it, we begin first with 
the Euler method itself, which is the oldest and most 
straightforward method of analysis, but is also 
relatively crude and inaccurate. 
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If the increment, h, in the Taylor series 
(21) is taken to be much less than unity, the first 
two terms will provide reasonably good accuracy. 
The solution at the neighboring point then becomes 


Y = Y + M h . (22) 

n+1 n n 


into the differential equation in order to find the 

value of M , , . This value is, in turn, used to 
n+1 

calculate a new and more accurate value of Y 

n+1 

from equation (23). The process can be repeated, 

using the new value for Y , , until the desired 

n+1 

accuracy is attained [19]. 


The solution is expanded to increasingly remote 
points from the initial point by the process 


Y 0 = Y + M h 
n+2 n+1 n+1 


Y = Y 0 + M 0 h , etc. 
n+3 n+2 n+2 


The trouble with the above method lies in the 
need to determine the derivatives of each variable 
involved. For the complex set of equations derived 
in the DERIVATION OF EQUATIONS section of 
this report, this is no easy matter. However, 
the same general averaging effect can be obtained 
in the following way, which, although less accurate, 
is more straightforward and easier to apply to the 
problem at hand. 


This process is known as the Euler method. 

While the first few points can be made 

sufficiently accurate by picking h small enough, 

as the solution progresses, the error will become 

increasingly larger. This occurs because the 

differential, M , (the slope of the solution curve 
n 

at x ) is used to calculate the solution Y . at 
n n+i 

x . The calculation of Y . will therefore be 
n+1 n+1 

in error by an amount equal to the difference in the 
slopes at x^ +1 and x^, multiplied by the incre- 
ment. Further the total error in the solution 

Y , is the sum of the error in calculating Y 
n+1 n+1 

and the errors in calculating all previous solutions. 

The fallacy in this method is discussed very well 

by Scarborough [18]. 

By using the modified Euler method, the 

cumulative error in the previous method is avoided. 

Here an average value of the differential is used in 

calculating each succeeding point. The extrapolated 

solution Y . is therefore of the form 
n+1 

Y =Y+-^/'M+M , \ h (23) 

n+1 n 2 ^ n n+1 / 


where M is the derivative of Y at x and 
n n 

M , is the derivative of Y at x . . This method 
n+1 n+1 

requires the use of a less accurate solution of 

Y such as equation (22), which is inserted 


First, we observe that the solution Y , can 

n-i 

be obtained by replacing h by -h wherever it 
appears in equation (21). Subtracting the resulting- 
equation from equation (21) , we obtain 


Y 


n+1 


Y +2Mh + -M”h J +... 
n-1 n 6 n 


If h is again taken to be less than unity and terms 
of order h 3 and higher are dropped, we obtain 


n+1 


Y + 2M h 
n-1 n 


(24) 


The accuracy of this formulation is much improved 
over that of the unmodified Euler method shown in 
equation (22). This is evident by the fact that the 
Euler method neglects all terms in the Taylor 
series of order h 2 or higher, whereas this method 
only drops terms of h 3 or higher. 

The extension of the above method to partial 
differential equations is rather straightforward, 
although stability is a much greater problem here, 
and some care must be exercised to ensure that the 
solution behaves in a stable manner [20]. 

The partial derivative of Y with respect to 

x at x and t, has the form 
n k 


dY \ __ Y n+l,k Y n-l,k 

Ox / 2Ax 

n 


(25) 


10 



which is of the same form as equation (24) . It 
would at first appear that the time derivative of Y 
could be written in the same form. However, this 
would require the storage of three arrays in the 
computer. Instead, the following form has proven 
to be more useful: 


/ dY \ _ Y n,k+i 2 Y n+1 , k + Y n-1, k 

\ at ; - At 

(26) 

Here, the second value in the finite difference has 
been averaged over two increments within the same 
time frame. This formulation will be stable pro- 
vided we choose At according to the relation 


At Ax/u 


(27) 


where u is the characteristic velocity involved. 
This is equivalent to requiring At to be smaller 
than the time required for a quantity to cross the 
Ax grid 121 J . 

For the solution at a boundary where n = 0 
and, therefore, the subscript (n-i) has no 
meaning, we require the less accurate forms. 


that the modified Euler method cannot be used to 
find the first two points of the solution. Two sets 
of difference formulas will, therefore, be required 
for any given equation. The first set will be based 
upon the relations (28) and (29) and will be used 
to find only the first two points in x for each time 
frame. The solution at all other points will be 
found from the second set of equations, which will 
be based on relations (25) and (26). 

First, consider the steady-state equations 
given in the Steady-State Simplification section. 

The form of the potential is given by equation (19) 
which we repeat below: 


d 2 P(.v) 

dy 2 


e P(y) 1 

n/ 1 - 2 a P (y) 


(19) 


For this equation, the accuracy of the unmodified 
Euler method given in equation (22) was found to 
be sufficient. An extension of equation (22) to 
the second derivative yields 


l2 Y - 2Y + Y 
d y n+1 n n-1 

dx 2 Ax 2 


(30) 


Applying equation (30) to equation (19) we find 


/ 3Y \ _ Y n,k+1 Y n,k 

\ 8t /n = ' At 


and 


/ QY \ Y n+ 1 , k Y n, k 

V )x / n AX 


( 28 ) 


(29) 


which are based on the Lin modified Euler method and 
will, in general, be used only for the first two terms 
of the solution. 


Finite Difference Equations 


P - 2P + p 
_ n+ 1 n n-1 

Ay 2 


n -l\/l-2aP 


solving Tor P . we obtain 
n+1 

P n+1 = 2P n ' P n-1 + ***** “ AyWl ~ ’ 

(31) 

which will provide a solution for the potential at all 
spatial points, with the exception of the first two 
at the origin. The ion velocity can be found from 
equation (16): 


The numerical equations can now be formulated 
by applying the appropriate relations in the foregoing 
section to the equations developed in the DERIVATION 
OF EQUATIONS section. Before going on, however, 
it should be pointed out, as in the Modified Euler 
Method section and as is apparent from equation (25) , 


P(y) 


1 

2 a 



(y) - 1 




(16) 


Solving this equation for V and writing the variables 

as functions of y , we obtain, 
n 
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V = n/ 1 - 2 a P . (32) 

n n 

From equation (18) we have the following relation 
for density: 


D (y) = i/V (y) (33) 

n n 

Finally, the steady-state potential on the plate, 
needed to give the solution the correct magnitude, 
can be found from equation (20) , which written as 
a function of y, is of the form 


D , - D , / V - V 

n, k+1 n, k n+l,k n,k 

AT n,k\ Ay 


_ v ( P n+l,k~ D n,k\ = 
n,k \ Ay / 


( 35 ) 


V - V 

n,k+l n, k 

AT 


/ V . , - V , \ 
y ( n+l,k n,kj 

” V n V Ay / 


+ a fc . =0 

n, k 


(36) 


P 0 = -J ln (D 0 Vq/V) 


(34) 


The equations (31), (32), (33), and (34) can 
now be programmed for a numerical solution to the 
steady-state case. This program is given in 
Appendix A. 

We consider now the time-dependent equations 
given in the Nondimensional Time-Dependent 
Equations section of this report. For convenience 
we repeat equations (11), (12), (13), and (14) 
below: 


jjgiy ill D(v T) lY.(y mII 
at ' ay 


- V(y,T) 


9D(y,T) 

3y 


0 


(11) 


A A 

E w , - E r 
n+1 k n, k 


V A 

^ E n,k Ay 


Ay 


D - e 
n, k 


n =y 


(37) 


and 


A A 

E , 4 - E , 

n,k+l n,k 

AT 


= D .V . - V e 
n, k n,k 


E „,k Ay 

n=y 


(38) 


Solving equations (35), (36), (37), and (38) for the 
appropriate variables, we find the following set of 
equations: 


D , = D , + 

n,k+l n,k 


n, k 


V - V 

n-f-l,k n, k 

Ay 


a - g y t : ' T - - V(y, T) 3 - V 3y ?y) + <* E(y, T) = 0 

( 12 ) 


A I E(y,T) dy 

8E ^ T) = D(y, T) - e y ' (13) 


+ V n,kV— ’O '/! AT < 3 ^ 


v = V ,+V I ■ n+1 ’ k EiiL ) AX 

n, k+1 n, k n,k \ Ay 


A 2 / E(y,T) dy 

— ( 8 ° x - j ' = D(y, T) V(y. T) - V e ° 


(14) 


Expressing these equations in the difference forms of 
equations (28) and (29) we obtain 


- a E AT 
n,k 


(40) 





oo 


n=y 


AE „, t « 


Ay 

(41) 
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AT 


A A 

E 0,k+1 = E 0,k 


■ - T ~ n=0 

*' D o,k v o,k- Ve 


*1 VkA 


AT 


/ 


(42) 


The above set of equations is used to obtain solutions 
to the variables at the plate boundary in all time 
frames. 

By applying the difference formulations in 
equations (25) and (26) to equations (11), (12), 
and (13), we obtain 


n, k+1 


: — ( v + V ^ 

2 y n+l,k n-1, k/ 


( V - V 

n+1, k n-1, k 

2Ay 

- of E .AT 
n, k 


( 47 ) 


A A 

E 4 , = E , 

n+l,k n-1 





(48) 


D , 4 - D . 

__ n ,k+l n,k 

AT 


- D 


V , , - V . . 1 
n+ 1 , k n-1, k 


n, k 


2Ay 


v ( ° n+1 ’ k " P n--1'A ) = o 

V n, k \ 2Ay / 


(43) 


V - V /V - V 

_ n, k+1 n,k I _ n +l,k n-l,k 

AT n,k \ 2Ay 


+ a E - 0 
n, k 


(44) 


A A 

E - E 

n+l,k n-l,k ^ n=y 

" — = D . - e 

2Ay n, k 




(45) 


where any quantity W , is defined as 
n, k 

1/2 ^ ^ ^ . Solving these equations 

for the appropriate variable we obtain the following: 


D n,k+1 2 [ °n+ 1 , k + °n-l,k 


V V 

D t _ n + 1 ? k - n~l_, k 
n,k \ 2Ay 


n, k 


D 4 , - D 4 , 

+ v l - n+1 > k 


2Ay 


AT 


( 46 ) 


This set of equations can be used in any time frame 
to find a solution to the problem at all points in 
space with the exception of those points on the plate. 

The above sets of equations, (39), (40), (41), 
and (42) and (46), (47), and (48), can now be 
programmed for a numerical solution of the time- 
dependent problem. This program is included in 
Appendix B. 


SOLUTION OF CASE I 

Application of Equations and 
Boundary Conditions 

In this first case to be treated it is assumed 
that the conducting plate initially has no electric 
potential. The plate is then allowed to float, its 
potential being determined by the ratio of the 
impinging ion and electron fluxes. The physical 
parameters governing this particular case are listed 
in Table 2. The boundary conditions imposed on the 
equations for this case are as follows: 

a. Initially, the parameters of the undisturbed 
plasma will apply at all points in space. On this 
basis, the electric field is zero everywhere while 
the nondimensional ion density and velocity are both 
unity. 

b. The ion density and velocity must at any time 
go to their undisturbed values (unity) and the electric 
field to zero at very large distances from the plate. 

c. The electric potential on the plate is given by 
the relation in equation (34) for the steady- state 
equations, and the electric field at the plate is given 
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TABLE 2. CASE I PARAMETERS 


Symbol 

Parameter 

Value 

- 

a 

Altitude (km) 

200 

N 

Ion-Electron Density 3 " (m“ 3 ) 

(3-50) x 10 10 

T 

e 

Temperature 3 " (K° ) 

450-800 

m. 

l 

Average Ion Mass 3 " (Amu) 

24 


b 1 


u 

Average Electron Velocity (km/ sec) 

43. 928 

V 

s 

Satellite Velocity 3 " (km/sec) 

7.780 

L 

Debye Length* 3 (cm) 

0.276 

1 

Potential Factor* 3 (volt/m) 

24. 973 

T 

Time Factor* 3 (sec) 

0.3548 x 10" 6 

. 


a. From Table 1 

b. Calculated Values 


at any time by equation (42) for the time- 
dependent equations. 

To determine the other parameters, such as 
non dimensional electron velocity, Debye length, etc., 
we must refer to specific values of the plasma 
parameters given in Table 1. With this information, 
the values_of the nondimen sional average electron 
velocity (V), Debye length (L), characteristic 
time (t), non dimen sional potential factor (£), 
and a factor that determines the effectiveness of 
the electric body force on the ions (n) , can be 
determined from relations given in the Nondimen- 
sional Time-Dependent Equations section. These 
values are given in Table 2 along with the data from 
which they were calculated. 

The steady -state program (Appendix A) is 
relatively simple, and, in principle, it can be made 
as accurate as desired by properly choosing the 
increment size, allowable error, and the magnitude 
of the assumed potential. Accuracy of more than 
four significant figures is unnecessary here, how- 
ever, since the general behavior is important and 
higher accuracy would not visibly alter plots of the 
data. This is the more accurate of the two pro- 
grams and should certainly be considered to give 
more reliable steady-state values. 


The program for the time -dependent solution 
(Appendix B) is considerably more complex, and 
an analysis of its accuracy is not so straighforward. 
Since it solves the partial differential equations, 
it has finite differences in both time and position. 

As a result, the error, which would normally be 
associated with a finite difference treatment of 
ordinary differential equations of one variable, is 
compounded by the second difference relation. 

This is further complicated by the restriction 
equation (27) places on the ratio of the increments 
of the two independent variables in order to attain 
a stable solution. This ratio was found to be 
extremely critical. It should also be noted that, 
because finite differences are used, the solution 
consists of small step-like increments rather than 
a continuous curve. The consequence of this can 
best be seen by taking the electric field as an 
example. Once the zero boundary condition far from 
the plate has been approached to within a critical 
value, the next solution step, rather than going 
smoothly to zero, crosses over into the positive 
region. Zero is missed because it is located only 
a fraction of an increment from the last step in the 
solution. However, once the electric field becomes 
positive, the solution becomes unstable and grows 
at an increasing rate. 
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The above effect is strictly a result of using 
finite difference approximations and has no physical 
foundation. However, to prevent its occurrence, 
it was necessary to force the electric field to go to 
zero whenever this situation took place. While this 
may seem arbitrary, it can be argued that this 
truncation does not affect the solution appreciably 
since the field is very close to zero whenever the 
instability occurs. 

Case I Results 

In the discussion to follow, the term steady- 
state will refer to the solution of the steady-state 
equations of the Steady-State Simplification section, 
while the solution of the time -dependent equations 
of the Nondimensional Time -Dependent Equations 
section will be described as asymptotically approach- 
ing equilibrium (hereafter referred to as quasi- 
equilibrium) . 

The plots of the quasi-equilibrium and steady- 
state electric potential given in Figure 2 are seen to 
be in very good agreement. The data points given 
in this figure are the results of the numerical 
solution of steady-state equations (31), (32), (33), 
and (34) where the curves result from the numeri- 
cal, quasi-equilibrium solution of time-dependent 
equations (39), (40), (41), (42), (46), (47), and 
(48) . The quasi-equilibrium values occur at 
t = 1.42 ps and we use the standard definition 
P 0 /e for sheath thickness. Hence, the quasi- 
equilibrium sheath thickness at 200 km is found 
to be 0. 33 cm. 



NON 0 1 HEMS I OML DISTANCE FROM MUTE MM 

Figure 2. Spatial distribution of electric 
potential at 200 km. 


The response of the ions to the electric field 
can be seen in Figure 3, which shows the develop- 
ment in time of both the ion velocity and density, 
along with the steady-state solution. At t = 1.42 ps, 
the time dependent solutions (denoted by curves) 
have essentially reached their equilibrium values, 
and the discrepancy between the quasi-equilibrium 
curve and the steady-state solution (shown as data 
points) is attributable to round-off error and the 
relatively large increment size used. This will be 
discussed in more detail at the end of this section. 

From the steady-state continuity equation we 
have the requirement DV = 1 everywhere in space. 



Figure 3. Ion density and velocity distribution 
in space and time at 200 km. 


Therefore, as a check on the accuracy of the 
numerical solutions, the quasi-equilibrium flux 
and steady-state flux at the plate (y = 0) and at 
the edge of the sheath (y = 1.2) were compared 
with unity. The quasi-equilibrium and steady- state 
flux differed from unity by 0, 000579 and 0. 0000152 
at y - 0 and 0. 000505 and 0. 0000520 at y = 1. 2 
respectively, thus showing the solutions to have good 
accuracy, with the steady-state solution being 
slightly better as expected. 

The growth of electric charge on the plate, 
given in Figure 4, indicates the charge behavior 
relatively early in the interaction, which primarily 
is due to electrons at this point. Initially, all 
electrons impinge on the plate and deposit charge 
there, but as the plate potential becomes more 
negative, increasing numbers of electrons are 
repelled and the rate of charge deposition decreases. 
The influence of the ions cannot be seen in the early 
stages of the interaction shown here. However, 
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Figure 4. Growth of surface charge on 
plate at 200 km. 

they have a very pronouned effect that takes place 
much later, which will be discussed in the following 
chapter. 

Figure 5 shows the development of the electric 
potential in time and space and illustrates very well 
the approach of the solution to an equilibrium value. 

An instability that occurred very early in the 
numerical solution of the time-dependent electric 
field is shown in Figures 6 and 7. This effect is 
strictly the result of error in the numerical cal- 
culations as is indicated by the difference in the 
two curves affected by changing increment sizes. 
Notice that decreasing the increment size (Figure 7) 
confines the instability to a smaller region in time 
and decreases its amplitude. In both cases the 
instability dies out in a small fraction of the 
relaxation time and is insignificant in the sheath 
region at all times. 

In Figure 8, the constant density contours in 
space and time give some idea of the lag in the 
response of the ions. The dashed line represents 
the position of the 1/e value (sheath) in space 
and time. 

Some additional comments should be made 
about the error involved in the time-dependent 
calculations. This is a result of both the numerical 
method and the computing machine used. Recall 
from the APPLICATION OF NUMERICAL METHODS 
section that the modified Euler method used herein 
is an approximation of the Taylor series where 
terms of h 3 or higher are dropped. Hence, to 
have good accuracy, the increment, h, must be 
made much less than unity. Therefore, in the 
present case, AT and Ay must be small, but, 
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Figure 5. Growth of electric potential 
in space at 200 km. 



Figure 6. Time dependence of the electric field 
for Ay = 0. 1 and AT - 0. 01 at 200 km. 
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Figure 7. Time dependence of the electric field for 
Ay = 0. 01 and AT = 0. 001 at 200 km. 
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Figure 8. Constant ion density contours 
at 200 km. 


in addition, we must satisfy the criteria for stability 
given in equation (27). Therefore, for an accurate, 
stable solution, we must have AT < Ay « 1, 

The above conditions, unfortunately, could not 
be completely complied with on the particular com- 
puting machine used. Smaller increments require 
a greater number of iterations, which in turn 
require more computer time. It was therefore 
necessary to use the rather large increments, 

Ay = 0.1 and AT = 0.01. 


The effect of using these increments can be 
seen if we consider the order of the terms dropped, 
h 3 , and the number of iterations. The missing 
terms are on the order of 0. 001 in space and 
0.000001 in time. The number of iterations 
required in space and time are 100 and 400, 
respectively. Hence, from this source alone, 
the accumulated error at equilibrium is on the 
order of 0. 1 and 0. 0004 for the two dimensions. 
This, coupled with machine roundoff error (which, 
using standard error analysis, has an accumulated 
value on the order of 0.0001) can easily account 
for the discrepancies in the quasi-equilibrium and 
steady-state solutions. 


SOLUTION OF CASE II 

Application of Equations and 
Boundary Conditions 

This second case treats the problem of a plate 
that is initially in equilibrium with the streaming- 
plasma for a given set of physical parameters. The 
plasma parameters then undergo a step change. 

The time-dependent equations derived in the 
DERIVATION OF EQUATIONS section are invoked 
to investigate the nature of the interaction that 
follows as the plate adjusts to the new plasma 
environment. 


This problem is described by the same set of 
time-dependent equations as Case I; namely 
equations (11), (12), and (13). For convenience, 
the physical constants corresponding to an altitude 
of 200 km and the steady -state distribution of the 
plasma parameters found in Case I are used as the 
boundary condition for this case at t = 0. After 
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the step change, the physical constants and plasma 
parameters corresponding to an altitude of 3000 km 
are used (Table 3) . We still require that at an 

TABLE 3. CASE U PARAMETERS 


Symbol 

Parameter 

Value 

- 

Altitude a (km) 

3000 

N 

Ion-Electron Density 3- (m -3 ) 

7 x 10 9 

T 

e 

Temperature a (K°) 

4000 

m. 

i 

Average Ion Mass 3- (Amu) 

4 

_ 

b 


u 

Average Electron Velocity (km/ sec) 

96.2 

V 

s 

Satellite Velocity 21 (km/ sec) 

6. 520 

L 

Debye Length* 3 (cm) 

5.2166 

4 

Potential Factor* 3 (volt/m) 

6. 6073 

T 

Time Factor* 3 (sec) 

8. 00 x i0" 6 


a. From Table 1 

b. Calculated Values 


infinite distance from the plate, the electric field 
vanish, the nondimensional ion velocity and density 
be unity, and at the surface of the plate, the electric 
field be given by equation (14) . 

The steady-state equations can be used as in 
Case I to obtain the equilibrium values of the first 
plasma state directly. This will provide the initial 
boundary conditions at t = 0 in a form that can be 
fed directly into the program for the time-dependent 
equations. 

The program for the solution of this case, which 
is a combination of the steady- state and time- 
dependent programs used in the previous chapter, 
is provided in Appendix C. Its operation is 
described sufficiently there. 


Case 1 1 Results 

The curve and data points given in Figure 9 
represent the quasi-equilibrium (t = 30, 5 ps) and 
steady-state solutions for the electric potential 
respectively. A comparison of these results 
shows good agreement between the time-dependent 



Figure 9, Spatial distribution of electric 
potential at 3000 km. 


and steady-state solutions for this case. Note that 
the electric potential has a greater magnitude than 
in Case I. Further, we find the sheath thickness 
to be 7.3 cm, which is also greater than in Case I. 
Therefore the effect of increasing the satellite 
altitude from 200 to 3000 km has apparently been 
to magnify the sheath dimension. 

The time and spatial distribution of the ion 
density and velocity are given in Figure 10. The 
quasi-equilibrium solutions are represented by 
curves and the steady- state solutions by data 
points. These parameters show much better agree- 
ment than in the previous case (Figure 3). This 
is probably because the transition in the plasma 
parameters is less drastic. Note that the change 
in ion density and velocity near the plate is about 
50 times greater than in Figure 3. This would be 
expected since we have already observed the electric 
field and electric potential to be considerably 
greater in this case. 

The effect of varying increment size on the 
spatial behavior of ion density and velocity is shown 
in Figure 11 for T - 0.35. Figure 12 shows the 
same effect in time for the ion density and velocity 
at the plate (y = 0) . Notice that in Figures 3 and 
10, the quasi-equilibrium solution decreases at a 
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Figure 10. Ion density and velocity spatial 
distribution at 3000 km. 



NON D 1 MENS I ONAL TIME (T) 


Figure 12. Effect of increment size on 
accuracy of time dependence. 
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Figure 11. Effect of increment size on 
accuracy of spatial distribution. 


slower rate than the steady- state solution. From 
Figure 11 it is apparent that one effect of reducing 
increment size is to increase the rate at which the 
time-dependent solutions decrease in space, thus 
tending to compensate for this discrepancy. The 
data shown in Figure 11 occurs rather early in the 
interaction (T - 0.35) so that the end effect on the 
qua si- equilibrium solution (T 3. 81) must be 
extrapolated. However, Figure 12 indicates that 
the effect of the decreased increment size will be 
more pronounced as the interaction progresses. 

The time dependence of the surface charge 
density on the plate (Fig. 13) indicates an 
interesting effect not observed in Case I. The 
surface charge density does not increase uniformly 
toward equilibrium as expected, but begins to 
decrease after an overshoot early in the interaction. 
The overshoot of surface charge density, which 
reaches its peak value at T = 0. 8, is the result of 
the electron interaction. The surface charge then 
begins to decrease slowly from the peak value as 
some of the negative charge on the plate is 
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Figure 13. Growth of surface charge on plate 
at 3000 km. 



ndndiuensional time in 


Figure 15. Growth of electric field at y = 4. 1 
and y = 5. 1 at 3000 km. 


effectively neutralized by the ions. The inertia 
of the ions is well demonstrated by this effect. 


Figures 14 and 15 give the time dependence of 
the electric field very early in the interaction at 
various spatial positions. The fluctuations at large 
distances from the plate that were observed in 
Case I (Fig. 6) are not found here; however, 
there is a slight overshoot at y = 4. 1 and y = 5. 1. 
Probably, whether merely an overshoot or fluctuations 
are observed depends upon the magnitude of the step 
functions used to represent changes in plasma param- 
eters. The gradual decrease late in the interaction 
is caused by ions, as previously discussed. In 
Figure 14 we see, however, that the point in time 
at which the electric field achieves its maximum 
value varies with distance from the plate. Apparently, 
the ions nearest the plate, which are consequently 



Figure 14. Growth of electric field at y = 0. 5 
and y = 1. 0 at 3000 km. 



DISTANCE FROM PLATE (y ) 


Figure 16. Growth of electric potential in space 
for transition from 200- to 3 000 -km orbit. 

exposed to the greatest electric potential, experience 
an acceleration earlier than ions located farther 
away. Hence the effect of the ions is felt later at 
increasing distance from the plate. 

Figure 16 gives the spatial dependence of the 
electric potential at various stages of the interaction, 
thus showing the approach of the solution to 
equilibrium. 
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DISCUSSION 


Conclusions 


Summary 

To gain some insight into the behavior of the 
ionospheric plasma directly ahead of the stagnation 
point on a satellite, a kinetic treatment of the prob- 
lem was carried out, based on two uncoupled, 
collisionless Boltzmann equations, corresponding 
to the two ionospheric species considered (ions 
and electrons) . The formulation of this treatment 
of the problem was based on the assumption of an 
infinite, flat, conducting plate representation of the 
satellite; a uniform ion velocity distribution in time 
and space; a Maxwellian distribution of electron 
velocity; and an instantaneous relaxation of the 
electrons. With these assumptions, the equations 
of transport derived from the two Boltzmann 
equations and a form of Gauss* law were used to 
obtain a set of four coupled differential equations 
which can be solved numerically to obtain the 
time-dependent plasma parameters. 

The first case to which these equations were 
applied consists of a grounded, conducting plate 
oriented perpendicularly to the plasma flow, which 
is allowed to float electrically at time t 0 . The 
reaction resulting as the plate moves toward an 
equilibrium state with the plasma is very similar 
to the type of reaction that occurs on board satellites 
when the potential at some point on the surface is 
periodically pulsed off and on by instrumentation, 
such as an emissive plasma probe. 


The second case, for which a solution was 
obtained, treats the problem in which the plate is 
initially taken to be at some preestablished 
equilibrium state. At time to, the plasma 
parameters undergo a step change and a time- 
dependent interaction follows as the plate moves 
toward a new equilibrium state compatible with 
its new plasma environment. 

This case is particularly applicable to the 
behavior of the sheath region ahead of a satellite 
as it passes over the terminator; through a radiation 
belt or anomaly; or, within certain limits, undergoes 
a change in orbital attitude [22] . 


The results of the two cases discussed above 
enable us to state the major events and parameters 
of the time -depen dent interaction and to draw the 
following conclusions: 

a. The plasma sheath thickness changes from 
0. 33 cm at 200 km to 7. 3 cm at 3000 km. At 

200 km the sheath grows from zero when the plate 
is grounded to its equilibrium value in 1.42 ps, 
whereas the transition from equilibrium at 200 km 
to an equilibrium state at 3000 km requires 30. 5 /^s. 
This would indicate that the relaxation time is 
inversely proportional to the plasma flux to the 
plate (both satellite velocity and plasma density 
decrease in the transition from 200 km to 3000 km) . 

b. When plasma parameters were varied step- 
wise, an initial overshoot was observed in the re- 
sponse of the plasma sheath and surface change on the 
plate. This overreaction of the plasma occurred very 
early and dampened out quickly, permitting the plate 
to continue toward equilibrium in a steady manner. 

c. One further event, which occurred before 
equilibrium was attained, is worth noting. The 
surface charge on the plate reached a peak negative 
value, G. 4 ps, after the step transition from param- 
eters at 200 km to those at 3000 km initiated the 
interaction. From this point, the surface charge 
moved toward an equilibrium value of less magnitude. 
This event illustrates the delay in the response of 
the ions and their influence on the interaction. 

d. The equilibrium potential on the plate was 
found to be 0. 06 V at 200 km and 0.46 V at 3000 km. 
This is seen to be in good agreement in both form 

(4 L - kT^/q) and magnitude with the predictions 

of YA. L. Al'Pert et al. [ 23 J . 

From the preceding discussion, it can be 
concluded that the general behavior of the plasma 
sheath directly ahead of the frontal stagnation point 
on a satellite is very close to that expected from 
physical arguments. In addition to confirming this 
general, qualitative behavior, the study presented 
here produces a more exact, quantitative description 
of the interaction that can be applied to many prob- 
lems with interesting and physically meaningful 
applications. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, October 9, 1969 
124-09-00-00-62 
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APPENDIX A 


A PROGRAM FOR THE STEADY-STATE SOLUTIONS 


This program is written specifically for the 
steady-state equations of the Steady-State Simplifi- 
cation section in the text. To facilitate a clear 
explanation of its operation, it will be divided into 
six parts. 

The physical parameters and constants in the 
first portion of the program will determine the exact 
physical problem to be solved. Therefore, the 
physical plasma parameters, temperature (XT), 
average ion mass (XMAS), and satellite velocity 
(XU) , must be properly assigned. The remaining 
parameters will be calculated from these, and the 
physical constants do not change. 

The second portion consists of only two state- 
ments. The first calculates the electric potential 
on the plate, and the second determines how closely 
this value must be approached. 

In the third part of the program, a small 
potential, VX1, is assumed to exist far from the 


plate. This potential, VXi, is defined by the first 
statement. The following statement defines the 
slope at this point, and the last statement calculates 
the potential at the next adjacent point. 

In the fourth section, the potential is calculated 
at points increasingly near the plate. When a cal- 
culated value exceeds VXO, the process is terminated. 
If the difference between the solution and VXO is 
DELT or less, the program goes to the next section; 
if not, VXI is changed slightly, and the process is 
repeated until the proper value of the potential at 
the plate is achieved. 

The next section merely reverses the potential 
and increment arrays so that the value at the plate 
becomes the first component. 

The last section calculates the values of ion 
density (XN) and velocity (W) at each point and 
prints out the results. 
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STEADY-STATE PROGRAM 


i 


Nr 3 
NMz ]. 

C PHYSICAL PARAMETERS AND CONSTANTS 

C ***************************************************** 

XTrlOOD. 

XMASr28,M 1 . 67 2 5 2* 10 . * * < -2 7) > 

XUr8. 0450*10. **3 
EMr9. 1091*10. 

XE r 1 *602*10* ** (—19) 

XK=1.38*in.**{-23> 

XYZ= «XK*XT >/ U 2* )* (3.1 41 59)*EM ) 

UE r S OR T ( XYZ) 

Ur (UE/XU) 

EPSXrXK*XT /( XM AS * ( XU * * 2 ) ) 

C ***************************************************** 

C VXD IS THE POTENTIAL ON THE PLATE 

C ***************************************************** 

VX Or ( . 5 ) * A LO G ( l./U) 

DELTr. OGI 

C ***************************************************** 

C VX1 IS SMALL POTENTIAL ASSUME AT SOME POINT IN SPACE 

C ***************************************************** 

vx ir-. ooon 1 

SLOPXrSQRT (1 »+E D SX)*VX 1 
DYr. 1 

YINC < 1 )rDY 
YINC (2 >rDY 
5 CONTINUE 
V( 1 ) r V X 1 

V(2)r5L0PX*DY+V( 1) 

C ***************************************************** 

C THE FOLLOWING PROGRAM EXPANDS THE SOLUTION RACK tqwARD 

C THE PLATE UNTIL VXD IS REACHFO. 

Z ***************************************************** 

Ir 1 

IT CONTINUE 
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YINC (1+2 > - D Y 

V( 1+2) = 2 .*V( 1+ 1) +( 0Y**2) *EXP (V (1+1 ) ) -C DY* + 2) 

1/SQRTC 1 • —2 **EPSX*V ( I + 1 )) —V ( I ) 

01 FX - A BS (VXO-V CI+2 ) ) 

IF (0 IF X-DFLT > 30*30. 15 

15 CONTINUE 

DI FX - V XO — V d+2 ) 

IF ( D IF X ) 16.1S.2C 

16 CONTINUE 
I- 1+ 1 

GO TO 10 
20 CONTINUE 

VX1 = VX1- .000 005 
GO TO 5 

c ***************************************************** 

C ARRAY IS REVERSED TO MAKE V(l)rVX n 

C ***************************************************** 

30 MMd+2 

MIDz (MM+ 1) /2 
DO 40 Jrl.MID 
K=MM+1 -J 
TTrYlNC( J) 

YINC (J >-YINC (K ) 

YINC (K )zTT 
T-V( J) 

V( J) -V (K ) 

40 V ( K ) -T 
M = MM 

C ***************************************************** 

C ION DENSITY AND VELOCITY A RF CALCULATED 

C ***************************************************** 

DO bH Iz 1 . M 

XN (I )zl. /SQRT< l.-2.*F°SX*V(I )) 

60 W(I)rl./XN(I> 

WRITE ( N.7C ) (I t YINC( I) *V (I ) t XN (T >♦ W( I) .1=1 .M .AIM) 

70 FORMAT (SX* INCREMENT POTENTIAL R EL . DENSITY 

1 R EL * VELOCITY V // / ( T6 tF<*. 5. 3£ 12 .7 )) 

STOP 

END 
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APPENDIX B 


A PROGRAM FOR THE CASE I TIME-DEPENDENT SOLUTION 


The program that follows is designed to solve 
the time-dependent equations of the Nondimensional 
Time-Dependent Equations section in the text. It 
is divided into six sections. 

The first section contains the various constants 
that determine step size and allowable error. Of 
particular interest are the constants DY and DT, 
which determine the space and time increment sizes 
respectively. Other constants of interest are: M, 
which determines how many values in the space 
curves are printed; MM, which determines the 
array size; DELP, DELD, and DELV, which are 
the allowable errors in the potential, ion density, and 
ion velocity solutions respectively; and LL, which 
determines how many of the time steps will be printed 
out. 

The second section consists of the appropriate 
values of the plasma parameters and physical con- 
stants. The values of the factors used to non- 
dimensionalize the variables are also calculated in 
this section and printed out. The names of these 
values are self-explanatory. 


The third section sets the initial conditions of 
the electric field, ion density, and velocity. These 
values form the first set of variables in time and are 
fed directly into the fourth section, which is set up 
to calculate the first two values of each time array. 
As pointed out in the APPLICATION OF NUMERICAL 
METHODS section, the first two values cannot be 
calculated by many of the more accurate numerical 
methods. 

Section five contains the bulk of the program. 
Here the remaining values of the array are calcu- 
lated and tested for convergence in space. The first 
value in the electric field array is then printed, and 
the solution moves to the next time increment. 

Section six contains the test for convergence in 
time and all the print statements for the variable 
arrays. Every LLth time step is printed out, plus 
whatever step the solution converges on. The con- 
vergence criteria in this section will not be satisfied 
unless E, D, and V have reached certain minimum 
differences in time. 
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CASE I. TIME -DEPENDENT PROGRAM 


DIMENSION EC 100*2) * D { 1 00 1 2 > * V C 1H V * 2 ) *P C 1 00 I 


t ****** * 


SECTION 1 


************************ 


*************** 


N- 6 

DYz.Ol 
DTr.001 
Mr 10 
MM r 1 00 

DEL-1. *10.** (-8) 

JlMrn 

DELP-.00001 
DIFPrDELP* 1. 

DEL D-* 000001 

DELVr. 000001 

KK r 1 

JJrlC 

J- J J“1 

LL-5 

Lro. 

ESlrO. 

Eso-r. 

W- 1 • 

***************************************************** 


SECTION 2 

***************************************************** 
XK z 1 .38054*10. **(-23 ) 

Tr 80 0 * 

EMz9. 1091*10. **(-31) 

XI Mr 24 .* (1 .57252*10. ** (-27 ) ) 

VS=7. 780*10. **(3) 

EVR-(XK*T) /( (2.)*(3. 14159) *£M) 

VArSGPT(EVR) 

AL PH A z (XK*T) /( XIM* (VS**2 ) ) 

Ur ( VA/ VS ) 

XNNz5. *1 0. **(11) 

E D Sz8. 8542*1 0. ** (-12 ) 

<3=1. 8021*10. **(-13) 

XY Zz EPS* XK*T/( (0**2) *XNN) 

XL r S OP T ( X Y 2 ) 




SQUIC-G*XNN*XL/EPS 
DTT-DT*XL/VS 
DX “D Y* XL 

WRITE! Nt 11 XLfSQUIGtPTT* DX »V A 

1 FORMAT (2X* 3HXL = t El 4. 5* 3X f6HSQUIGrf E14. 5t 3X *4HDTTr, El 4. 5# 3X .3HDX- t 
lE14.S»3Xt3HVAr,E14.5) 

***************************************************** 

SECTION 3 

***************************************************** 

DO ID I- 1 1 MM 
E( It 1)=0. 

Dt i* n n . 

ID V« It 1) =1 • 

***************************************************** 

SECTION 4 

***************************************************** 

2 0 EUt2)=Efl»l>+(Dlltll*VCltl)-U*EX*»C7.*ESl*DYn*0T 
2 ? CONTINUE 

D( 1 t 2) -D ( 1 f 1 ) + D( 1. I) ♦( V( 2f I) -V (1 f 1 ) ) * ( DT/OY) +V (1 f 1 (D (2 f 1 )-0t 1 * 1) 
1 ) * (DT/DY ) 

V( 1 .2) -V (1 • 1 )4*V( If 1) *( V( 2. 1) -V <1 f 1 n *CDT/OY) -ALt>HA*E (1 1 1 ) *DT 
E 1 -E ( 1 t2 ) * (*0001 ) 

EE -E ( 1 t 2 ) 

1-1 

IE (KK-2) 24*26,26 

24 CONTINUE 
APROXrO. 

GO TO 28 
26 CONTINUE 

APR0X-ES1+ (ES1-ESC) 

28 CONTINUE 

ES2--ABS ( APROX-tE) 

E( I+lf2)-E (I t2)+(0(I*2 )-EXP(E^2*PY) ) * D Y 
EE-EE+ ( E ( I +1 *21 — E(It 2) (T + l .5 ) 

J-I + 2 

0<It?)-(D(I— ltl)+D(I + l*l ))/?#+( n (Ttl )*(V(I+1 *1 )-VtI-Ifl) )+VtIfl) * 

1 (D ( 1 + 1 » 1 )-D( I- 1 1 1) )) *HT/ C ? .*0Y ) 

V(I»?)r(VU-lf l)+va+lfl l)/?. + V(If l)*{V(I+lf 1 )’V(I- 1 » 1 ))*DT/ 

1 (?.*DY)-AL PH A* CE ( 1-1 »1 J + E ( I + 1 * 1) ) * C T / 2 • 

3D CONTINUE 

**************************** **************** ** ******* 


SECTION 5 


***************************************************** 

ES2--ABS ( APROX-EE) 

E( 1+ It 2) -E <1-1 t2) + 2. *( DC It 2) -EXP C^S2*DY) MPY 
TE 5 T -E C I + 1 t 2 ) 

IF (TEST) 3 4 t 34 * 3 2 
32 CONTINUE 
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EC 1+1* 2) =0. 

34 CONTINUE 

E F r £ E + (E ( I +1 *2 ) — fE < I* 2 > >/?.+ECI+l»?) 

1 = 1+1 

TE ST -E C I +1 1 1 ) 

IF (TEST) 36* 35 *36 

35 CONTINUE 
DC 1.2) =1 . 

V< 1*2) =1 • 

GO TO 33 

36 CONTINUE 

D(I* 2 )=CD(I-l*l)+D(I*l*in/ 2 .+ (D«Ttl>*CV(I+l»l)-V(I-lfin + VCTtl)* 
1 CDCI+1 * 1 I-DCI-ltl) ))*QT/(2**OY) 

VCI.?>=(V(I-l*l>+V<I+l*l))/2. + V(I*l>*<VtI+l*l>-V(I-l*in*DT/ 

1 C2 .*0Y >- ALPHA* CE C I -l *1 ) + E ( 1+ 1 . 1 ) ) *DT/2 . 

38 CONTINUE 
NNrMM- 1 

IF(I-NN) 40*50*50 
4 0 CONTINUE 
GO TO 30 
50 CONTINUE 

EC 1+ It 21 = E (1-1 t2) + ?.*( D< It 2) -EX* (FS2*DY) 1*DY 
TESTrE <1+1 .2) 

IFCTEST) 5 6 t 56 *54 
54 CONTINUE 
ECI+It 2)=0. 

56 CONTINUE 

EE = E F + (ECT+1 *2)-E(I*2) )/2.+ECI+lt?l 
DC 1+ 1 * 2 ) =D < I *2 ) 

VCI+1*2)=V(I*2> 

MMrl+l 

DIFArABS CAPPOX-EE) 

DELTrABS CEE*. 01) 

IF CD IF A- OE LT ) 70*70.60 
oD CONTINUE 

DIFAr < APROX-EE ) 

APROX= APROX-OIFA/W 
W=W + • 5 * W 
1=1 

EE = E { 1 ,2 ) 

GO TC 23 

70 CONTINUE 
W= 1 • 5 

WRITF CN.8P) KK.EC1 .2) 

30 FORMAT (6Xt 22HELECTRIC FIELD AT TINE * 2X »I4» 2X *7 HE (1 *2 )-»£14 »5) 
***************************************************** 

SECTION 6 

♦♦♦♦a************************************************ 

L = L + 1 

IF(L-LL) 85*85.84 
34 CONTINUE 
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OIFP-ABS (EE— ESI) 

DIF0-ABS(0(1 *2 ) - D ( It I) ) 

OIFV = ABS (V Cl • 2 )-V( 1*1)) 

IF (DIFP-DELP ) 300*300.85 
?0G CONTINUE 

IF(DIFD-DELO ) 310. 310* 85 
310 CONTINUE 

IF (DIFV-DELV ) 90.90.85 
85 CONTINUE 
J- J+l 

IF(J-JJ) 160.90.90 
90 CONTINUE 
J-0 

WRITECN* 100) KK 

1?0 FORMAT C//2X. 25HPLASMA PARAMETERS AT TIME. 2X. 14) 

NN -MM- 1 
SUMrEE 

DO 110 I r 1 *NN 
P( I) - 5 UM *0 Y 

SUM- SU M- (F(ItZ> + ECI+l«2))/2. 

IF (SUM ) 110. 110* 104 

104 CONTINUE 
SUMrr. 

110 CONTINUE 

P ( MM ) - P ( NN ) 

WRITE( N* 120) (I.P(I).m*MM.M) 

no FORMAT (/ //2X *9HP0TENTI AL// (6 X* If. *E 16 .7 *16* El 6* 7* 16 »E 16 .7 .1 6. El 6. 7, 
1I6.E1S.7)) 

WRITE? N. 130) CI*ECI.2)»I=1»MM,M| 

130 FORMAT (///2X*1 4 HEL EC TR IC FI« r LD//(6X,I6,E16.7.I6.E16n,I&.F16.7.I6. 
1 E 1 6 • 7 * I6*F16*7) ) 

WRITE(N* 140) (I. 0(1.2) .Irl.MM.M) 

! 4 C FORMAT C///2X.18HI0N NUMBER DFNSTTY//(6X.I6*E16.7*I6*Elo.7,I6.F16.7 
1 *I6.Elo*7* 16 *E16*7 I ) 

WRITEtN* 150) (I*V(I*2).I=1.MM,M) 

I 5 C FORMAT (///2X .20HAVERAGE ION VELOCITY// (6X* 16 . E 16 *7 *I6»E16* 7. IS * 

1 E 1 6 • 7 * I6»F16.7*I6*E16«7) ) 



WRITE ( N. 155) 


155 

FORMAT (/// ) 



IF (L-LL) 160*160.156 

156 

CONTINUE 
IF (DIFP-DELP ) 

3 ? 0 » 32 0. 330 

320 

CONTINUE 
WRITE ( N* 31 5) 


315 

FORMAT ( 3 X . 11HP 
JI M- JI M + 1 

CONVERGES/) 

3 30 

CONTINUE 
IF (DIFD-DELD ) 

340, 340. 350 

340 

CONTINUE 
WRITE! N* 345) 


345 

FORMAT ( 3 X . 11HD 
JI M - JI M + 1 

CONVERGES/) 

3 5 0 

CONTINUE 
IF (DIF V-DELV ) 

360* 360* 160 
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?6 C CONTINUE 

WRITE(N?365) 

365 FORMAT ( 3 X » 11 HV CONVERGES/) 

JIMr JIM+1 

IF (JIM-3 ) 16 0? 190? 190 
150 CONTINUE 

DO 180 I r 1 • MM 
E(Ifl)-E(I»2) 

D(I?1)-0(I?2) 

130 V(l»l)-va.2) 

ES0-F.S1 
Esirre 
EE -Q • 

KKrKK+ 1 
JIM-0 
GO TO 20 
190 CONTINUE 

WRITE ( N? ?00) 

700 FORMAT ( 2 Xt 22HE OUILIBRI UM CONDITIONS) 
STOP 
END 



APPENDIX C 


A PROGRAM FOR THE CASE 1 1 TIME-DEPENDENT SOLUTION 


The program in this appendix is designed to 
solve the general time- dependent problem treated 
in case II. In this form, the program can solve the 
time-dependent equations developed in the 
DERIVATION OF EQUATIONS section for any of 
the various applications pointed out in the DIS- 
CUSSION. The program is divided into two sections 
for the following discussion. 

The first section essentially consists of a 
revised form of the program presented in Appendix 
A. Here the variable arrays, which were one- 
dimensional in Appendix A, are written in a two- 
dimensional form so that they can be applied 
directly to the time-dependent portion of the pro- 
gram in the form of initial conditions. In this first 
section, the physical constants and plasma param- 
eters constitute the initial conditions, or the initial 
plasma state. The constants and plasma parameters 
at the first of section two will constitute the final 
plasma state. 


The second section of this program is almost 
identical to the program presented in Appendix B. 
Therefore, we will point out and explain the 
differences in the two here. 

The form and purpose of the various constants 
are the same here as in Appendix A. However, the 
values of the plasma parameters have been changed 
in order to formulate the problem described in 
the SOLUTION OF CASE II section. 


The part of the program referred to as section 
two in Appendix B has been omitted here. This 
section provided the boundary conditions, which 
are calculated in section one of this program. 


The remainder of the program is identical to 
that presented in Appendix B and will therefore not 
be discussed further here. 



CASE n. TIME-DEPENDENT PROGRAM 


DIMENSION YINCCIOOni 

0 1 MENS ION E( 100C*2) t P< 1000 t? ) • V( ICTOt 2 1 f P( 10CC ) 


********************** ********* 


SECTION 1 


t ****** i 


c***** ****************** 


PLASMA STATE 1 PARAMETERS 

N-b 
Mr 10 
KK =0 

XU = 7. 78*10. **{ 3) 

XT=8CC . 

EMz9. 1091*10. **(-31) 

XE - 1 .802*10.** (—19) 

XK-1 .38*lC.**(-23> 

XMAS=24*(1. 67252*10. **(-27)) 

DELT=. 001 

EPSXrXK*XT/(XMAS*(XU**2) > 

D 0 - 1 . 

V0 = 1 . 

XY2-(XK*XT)/((2. )*(3.14159)* r M) 

UE=SGRT( XYZ) 

U= (UE/XU ) 

VXO = (.5) * A LO G ( 1. /U> 

VX 1=-. 00C05 
1039 CONTINUE 

SL0PX=SQRT(1.+EPSX)+VX1 
QY = . 1 

DO 1001 1=1*110 
P(I)=0. 

1T91 CONTINUE 
P( 1 ) =VX1 

R(2)=SL0PX*DY+P( 1) 

1=1 

lriq CONTINUE 

YINC(I+2)=QY 

P(I + 2)=2.*P(I+1)+(0Y**2) *EXP(P(I+1 ) ) - ( D Y * * 2 ) /SORT ( 1. -2 .*FPSX 
1 *P ( I +1 ) ) -P (I > 

DIFX=AB5 (VXO-P ( I+? ) ) 

IF(OIFX-DELT) 1C?0#1030»1C15 
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1015 CONTINUE 

DIF X~VXO —P (1+2 ) 

IF ( 0 IF X ) 1016# 1015 »1P2C 

1 D 1 6 CONTINUE 
1 = 1+1 

GO TO 1010 
1020 CONTINUE 

VX1=VX1-.0OCOC1 
GO TC 1005 
1030 MM-I+? 

MID=CMM+l>/2 
DO 1040 JrltMID 
K=MM+1 -J 
TT = Y INC C J) 

YINC ( J ) -YINC (K ) 

YINC CK ) = TT 
T - P ( J ) 

PC J ) zo (K ) 

1040 PCK)=T 

DO 1060 1=1# MM 

V<I#l)=S0RT(l.-2.*EPSX*PCin 
1050 DC I* 1) =1 ./VC I» 1) 

WRITECNtl2) KK 

12 FORMAT C//2X# 25H°LASM A PARAMETERS AT TIME»2X#I4) 

**+* 4 +*+*+++**+ + *+ + ***++***++*++***+***********+***+* 

SECTION 2 

***************************************************** 

PLASMA STATE 2 PARAMETERS 

DY=1C. 

DT = • CD 1 
MMzSPO 

DEL= 1. *1 0. ** C — 8 ) 
jlM = r 
DELP = .CC5 
DELD=. CGD001 
DELV=. OGOOOl 
L= 0 
LL = 5 
KK = 1 
JJ = 1^ 

J= JJ-1 
ES1=C. 

ES 0 = 0 • 

W= 1 • 

XK = 1 .38054*10. ** C - 2 3 ) 

T= 4 0 OG . 

EM = 9. 1091*10. **C-31) 



XI Mr 4. *{ 1.67252* 10.** C -27) ) 

VS =6 .520*10. ** (3) 

EVRr CXK*T) /( C2.) *C 3. 14 159) *EM> 
VA=S0RT(EVR) 

ALPHA=(XK*T)/(XIM*(VS**2>> 

U = i VA/VS ) 


XNN=7. *1 C. ** (9 ) 

EPS=8. 854?*1 C. ** (- 12 ) 

0=1. £021*10. **(-19) 

XYZ = EPS*XK*T/< (Q**2> *XNN > 

XLrSGPTC XY?) 

SQUI6=Q*XNN* XL/EPS 
DTT=GT*XL/ VS 
DXrOY*XL 

WRITE ( N t 1 > XL»SQUIG» OTT* DX 

1 FORMAT (2X*3HXL=*E14.5» 3X * 6 HS QU 19 = * E 1 4 . 5 1 3X *4 HD TT = t El 4 . 5 1 3X 1 3 MO X= » 
1E14.5) 

WRITE ( N» 16 ) (ItP(I)f 

IB FORMAT C///2X »9HP0TENTI AL// (6 X* 16 t^lS .7 *I6t El 6. 7t 16 tE16 .7 *I6» E16. 7* 
1 IS »£ 1 6 .7 ) ) 

1=1 

E< I* 1) =- (P (2 )~P( 1> ) / DY 
DO 10 I r ?• MM 

E( It 1) =- (P (I +1 )-P( I- 1) )/(2.*DY> 

10 CONTINUE 

WRITE ( Nt 17 ) (I tE Cl *1 >, ir 1, MM»M ) 

1 7 FORMAT C///2X »1 4H ELECTRIC F IFLD // CS X » I6?E16.7 tI6«E16.7t 16 »E 16 .7 *16* 
1E16.7# 16 *£16.7)1 
WRITE ( N* 18 ) Cl .0 (I .1 ) f 1= 1» MM f M j 

18 FORMAT C///2X rl8HI0N NUMBER DFNSTTY// (6X* 16 *£16 .7 # 16* El 6. 7» 16 *E16 .7 
1 tI6»E16.7t I6.E16.7) ) 

WRITECNf 19) (ItVCI.l )• I=lt MM. m) 

IS FORMAT C///2X *20HAVER AGE ION VE LO C I T Y / / ( 6 X * 16 *E 16 .7 1 1 6 t E 1 6 . 7t IS t 
IZlB.lt 1 6 tE 16. 7,16 » El 6. 7) > 

WRITEC Nt 1000 ) 

IfOC FORMAT (/// ) 

20 E( 1 » 2 ) =E (1 *1 ) + ( 0 ( 1 » 1 )*V( 1? 1) -U*FX D C2.*ES1*DY ))*0T 
GO TO 23 

21 CONTINUE 

E( It 2) =E Cl »1 ) 

23 CONTINUE 

DC 1 1 2) =D (1 rl ) + Dt It 1) * ( VC 2, 1) “V (1 1 1 ) ) *<DT/OY) +V Cl ,1 )* CD (2 1 1 >-0< It 1) 
1 )* (DT/OY ) 

V ( 1 » 2 ) ~ V ( 1 *1 ) * V ( 1*1) *(V(?t 1 ) -V (1 tl ) )*<OT/DY) -ALPHA*Etl»l )*DT 
£ 1 =E < 1 »2)*(.C0ni ) 

EE=E ( 1 * 2 ) 

1=1 

IFCKK-2) 24,26,26 

24 CONTINUE 
APROX=D. 
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GO TO 28 
26 CONTINUE 

APROX-ES 1+ (ESI ~E SO ) 

28 CONTINUE 

ES2=-ABS (APROX-EE) 

IF(V(1»2)-V( 1*1) ) 40 Ot 400*410 
4CC CONTINUE 

E( I+lt2)“E ( I + 1 *1 ) 

GO TO 420 
4 1 C CONTINUE 

E( I+lt 2) rE (I *2 ) + (0( I *2 ) — EX P ( FS 2* 0 Y ) ) *D Y 
4?0 CONTINUE 

EE=EE+ (E (I +1 *2)-E(I*2> >/?.+E (1 + 1*2) 

I - 1+ 1 

0(I#?)-(D(I-“l*l)+D(I+ltl )>/2.+(D(T*l )*(V (1+1*1 )—V(I-l»l) ) + V(X*l) 
I (0(1+1 *1 ) - D ( 1-1*1) )) *0 T/ ( 2 .* 0 Y ) 

V( 1*2) - (V( I- 1 * D+VCI + l »I )) /2«+V(T* 1)*(V(I+1* 1) -V (1-1*1 ) ) +0 T/ 
l(2.*0Y)-ALPHA*(E(T-lfl)+E(I+l*l))*DT/2. 

3 0 CONTINUE 

ES2“-A3S (APROX-EE) 

IF ( V ( 1 *2 )- V( 1* 1) ) 4 3 0* 43 Of 44 0 

430 CONTINUE 

E( 1+1* 2) =F (1+1 *1 ) 

GO TO 450 
44D CONTINUE 

E(I+I#2)=E(I-I.2) + 2**(D<If2)-EXP<ES2*0Y)>*DY 
450 CONTINUE 

TEST -E (1+1 *2 ) 

IF(TEST) 34*34*3? 

32 CONTINUE 

E ( I+It 2) =0. 

34 CONTINUE 

EE=EE+(E(I+lt2>-E(I*2> )/ 2. +E ( I +1 * 2 ) 

1=1+1 

TESTrE (1+1*1) 

IF(TEST) 36*35*36 

35 CONTINUE 
0( It 2) =1 . 

V< It 2) =1 • 

GO TC 33 

36 CONTINUE 

G(I*?)-(D(I— 1* 1) +0(1+1 *1 ))/?•+ (0(1*1 ) * ( V ( I +1 »1 )— V(I—1*11 ) + V ( 1 1 1) 
1 (0(1+1 *1 )— 0 t I- 1 * 1) ) ) *0T/ ( 2 •* OY ) 

V ( 1*2) = ( V ( I— It 1) + V ( I + 1 *1 )) /2 • + V( I* 1)*(V(I+1* 1) -V (I-l *1 ))*0T/ 

1 (2**DY)— ALPHA* ( E ( I - 1 tl )+E( I+lt 1) ) * D T / 2 • 

38 CONTINUE 
NN - M M- 1 

IF ( I — N N ) 4 0 1 50*50 



40 CONTINUE 
GO TO 30 
50 CONTINUE 

IF ( V 1 1 *2 )— V ( It 1 1 > 4 6 Ot 46 0* 47 n 
460 CONTINUE 

E(I+1.2)=E(I+I»I) 

GO TO 480 
47C CONTINUE 

EC I+lt 2) rF CI-1 .2 ) + 2. D( It 2) -CXP <FS?*DY) ) *0Y 
430 CONTINUE 

TESTrE <1+1 *2} 

IF (TEST) 56*56*54 
54 CONTINUE 

El 1+ 1 1 2 ) “C • 

56 CONTINUE 

EE-EF+ (E (1+1 *2)— E(I* 2) >/ 2. +E CT +1 t ’ 1 
0(1+ 1*2) “0(1*2) 

V( I+1*2)=V (I *2) 

MMrI+l 

DIFArASS (APROX-EE) 

DE L T - A BS (EE* .0 1 ) 

IF ( D IF A-PFLT ) 70 *70, 60 
6 C CONTINUE 

OIF A:( APROX-EE ) 

APROX-APPOX-DIFA/W 
W- W + .5 * W 
I - 1 

EE-E ( 1 *2 ) 

GO TO 28 

70 CONTINUE 
W= 1 • 5 

WRITE ( Nt 80 ) KK t E ( 1 *2 ) 

30 FORMAT (6Xt 22HELECTRIC FIELP AT TTwr,2X *14, ?X *7 HE (1 »2)-»EI4.5) 
L - L + 1 

IF(L-LL) 35*85*34 

84 CONTINUE 
QIFpzABS (EE-ES1) 

DIF0 = A8S in ( 1 »2 )-0( It I ) ) 

DIFV-ABS (V ( 1 *2 ) — VC 1*1)) 

IF(DIFP-DELP) 300*300*85 
7 0C CONTINUE 

IF (OIFD-DELO ) 31 0* 31 0* 85 
7 1 C CONTINUE 

IF (OIFV-DFLV ) 30*30*85 

85 CONTINUE 
J= J+ 1 

IF(J-JJ) 160*90*90 
9 G CONTINUE 
Jr 0 

WRITE(N* 100) KK 
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lHD FORMAT <//2Xt 25HPLAS«A PARAMETERS AT TIME. 2X, 14) 

NN=MM- 1 
SUMrFE 

DO 1 15 IrltNN 
P( I) -SUM*OY 

SUMrSUM- (E (I ,2 )+E( 1+ 1, 2) )/2. 

IF (SUM > llOt 110. 104 
104 CONTINUE 
SUMrO. 

lie CONTINUE 

PC MM ) = P{ NN ) 

WR I T E ( N . 120) ( ItP(I) il-l t M M f M) 

120 FORMAT C///2X t9HP0TEMTI AL// «6X. Ifi #F16.7t I6t El 6. 7. 16 ?E16 .7 . 16. E16. 7t 

IIS. £16. 7)) 

WR ITE(N» 130) ( If E(If 2) *1-1 tMMfM) 

130 FORMAT <///2X .14HELFCTRIC F IELO // < 6 X . 16 » E 16 .7 . 1 S t FI 6, 7 « 16 ,E 16 .7 1 1 6. 
1 E 1 6 • 7 1 I6.E16.71 ) 

WRITE ( Nt 14 0) ( I . 0 ( I» 2) *1-1 tMMt M) 

14G FORMAT <///2X *18HI0N NUMBER 0 EN S I T Y / / ( 6 X t 16 tE 16 .7 t T & . E 1 6 . 7. 16 f E 16 .7 
1 tl6tE16.7*I6»E16.7>) 

WRITE( Nt 150) (ItVIIt?) tl-l tMMfM) 

ISC. FORMAT (///2X #2GH AVER AGE ION VFL0CITY//(6X.I6.E16.7.r&.E16.7»IS. 
lE16«7»I6tE16.7»IG»E16.7) ) 

WRITE ( N« 155) 

155 FORMAT (/// ) 

IF (L-LL) 1 60 *1 60 . 1 56 
ISC CONTINUE 

CONVERSANCE TEST 

IF (OIF P-OEL p ) 3 2 C r 32 P * 33 P 

72C CONTINUE 

WRITF ( Nf 31 5) 

315 FORMAT (3X. 1 1 HP CONVERGES/) 

JTMr JIM+ 1 
330 CONTINUE 

IF(OIFD-OELO) 340*340. 35 0 
340 CONTINUE 

WRITE ( N v 7 4 5 ) 

345 FORMAT (3X» 11HD CONVERGES/) 

JIMz JIM+ 1 
350 CONTINUE 

IF ( Q I c V- DE L V ) 3&C. 3SP. 16H 

^0 CONTINUE 

WRITE ( N» 355) 

355 FORMAT (3X» 11HV CONVERGES/) 

JIM-JIM+ 1 

IF (JIM-3 ) 16CI. ISO. 1 5 C 

1 6 C CONTINUE 
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RE! SE T FOR NEXT TIME ITERATION 


Vl'Vd fi ) 

DO 1 8 C IU.MM 
t{I,l):£(I»?) 

D(If 1)=0(I«2) 

1?0 V(I» l)=V(It2) 

ESO=ESl 
ES 1 -EE 
EE -0 . 

KK ~K K+ 1 
J I M - 0 

IF C V ( 1 *2 > — V 1 ) 21 » 2 1 * 20 
190 CONTINUE 

WRITE ( N ♦ 20 0) 

200 FORMAT <2X* 22HEGUILI8RIUM CONDITIONS) 
STOP 
END 
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